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Abstract 


The increased flexibility of long endurance aircraft having high aspect ratio wings necessitates at- 
tention to gust response and perhaps the incorporation of gust load alleviation. The design of civil 
transport aircraft with a high aspect ratio strut or truss braced wing furthermore requires gust re- 
sponse analysis in the transonic cruise range. This requirement motivates the use of high fidelity 
nonlinear computational fluid dynamics (CFD) for gust response analysis. This paper presents the 
implementation of gust modeling capability in the CFD code FUN3D. The gust capability is verified 
by computing the response of an airfoil to a sharp edged gust. This result is compared with the theo- 
retical result. The present simulations will be compared with other CFD gust simulations. This paper 
also serves as a users manual for FUN3D gust analyses using a variety of gust profiles. Finally, the 
development of an auto-regressive moving-average (ARMA) reduced order gust model using a gust 
with a Gaussian profile in the FUN3D code is presented. ARMA simulated results of a sequence of 
one-minus-cosine gusts is shown to compare well with the same gust profile computed with FUN3D. 
Proper orthogonal decomposition (POD) is combined with the ARMA modeling technique to predict 
the time varying pressure coefficient increment distribution due to a novel gust profile. The aeroelas- 
tic response of a pitch/plunge airfoil to a gust environment is computed with a reduced order model, 
and compared with a direct simulation of the system in the FUN3D code. The two results arc found 
to agree very well. 
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Nomenclature 


a 

Cloo* 

A 

b 

C* 

c p 



/ * * r * r * 

x iJy iJz 
8 

G g _ 

i,J,k 
^8 cos 
Lg exp 
Lg omc 

Lgsin 

T * 

^ 8 cos 

J * 

^ 8 exp 

I * 

omc 
Lg* sin 

Lr* 

Mcc 

M 
bf hi 

N s 

N t 

NpOD 

p 

Q 


s 

t 

t* 

tr efcos 
t re f exp 
^ re f omc 
tre f sin 

IJ * 

Ly oo 

U gcos 
^8 exp 
^8 omc 


Autoregressive (AR) coefficients 
Free stream speed of sound 
State matrix 

Moving average (MA) coefficients 

Dimensional airfoil chord 

Surface pressure coefficient 

Surface pressure coefficient increment due to gust 

Roger approximation coefficient matrix 

Roger approximation coefficient matrix 

Arbitrary functions defining dimensinal x,y,z components of gust velocity 
Generalized displacement 
Generalized force due to gust 
Cartesian unit vectors 

Cosine gust wave length in grid units ( L gcos = L g * co jL R *) 

Gaussian gust length in grid units {L g = L g * exp /L R *) 

One-minus-cosine gust length in grid units (L gomc = L g * omc /L R *) 

Sine gust wave length in grid units (L gsin = L g * sin /L R ) 

Dimensional cosine gust wave length 
Dimensional Gaussian gust length 
Dimensional one-minus-cosine gust length 
Dimensional sine gust wave length 

Reference length for FUN3D spatial nondimensionalization 

Free stream Mach number 

Structural and aerodynamic mass matrix 

Number of structure modes 

Number of surface points 

Number of time steps 

Number of proper orthogonal decomposition modes 
Number of autoregressive (AR) terms 
Number of moving average (MA) terms is Q + 1 
Free stream dynamic pressure, psf 
Aerodynamic time, ( S = 2U„*t* /C*) 

Integration matrix 

Nondimensional time ( t = t*a 00 * /L R *) 

Dimensional time 

Cosine gust nondimensional reference time 
Gaussian gust nondimensional reference time 
One-minus-cosine gust nondimensional reference time 
Sine gust nondimensional reference time 
Dimensional free stream velocity 

Cosine gust nondimensional x-direction velocity magnitude 
Gaussian gust nondimensional x-direction velocity magnitude 
One-minus-cosine gust nondimensional x-direction velocity magnitude 
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U 8sin 
U 8 * 

^8 cos 
^8 exp 
^ Some 
^8 sin 

V 

W 8cos 
W Sexp 
^ Some 

W 8sin 

* 

w g 

*,y,z 

x*,y*,z* 

*T,yT,ZT 

a 

P 

P 

7 

A 

0 


X 

a 


Sine gust nondimensional ^-direction velocity magnitude 
Dimensional ^-direction gust velocity magnitude 
Cosine gust nondimensional y-direction velocity magnitude 
Gaussian gust nondimensional y-direction velocity magnitude 
One-minus-cosine gust nondimensional y-direction velocity magnitude 
Sine gust nondimensional y-direction velocity magnitude 
Dimensional y-direction gust velocity magnitude 
Cosine gust nondimensional z-direction velocity magnitude 
Gaussian gust nondimensional z-direction velocity magnitude 
One-minus-cosine gust nondimensional z-direction velocity magnitude 
Sine gust nondimensional z-direction velocity magnitude 
Dimensional z-direction gust velocity magnitude 

Location of point in nondimensional Cartesian coordinates, (e.g. x = x*/Lr*) 

Location of point in dimensional Cartesian coordinates 

Nondimensional grid speed metrics 

Modified nondimensional grid speed metrics 

Angle of attack, degrees 

Vector of coefficients in proper orthogonal decomposition analysis 

Predicted value of /3 

Roger approximation lag root matrix 

Structural and aerodynamic damping matrix 

Structure eigenvector matrix 

Vector of training values of dynamic component of surface 
nodal pressure coefficients 
Predicted value of S 
State variable matrix 

Structural and aerodynamic stiffness matrix 


1 Introduction 

Aerodynamic efficiency increase and drag reduction are key NASA Subsonic Fixed Wing Program 
goals. [1] A component of that research program investigates the truss braced wing (TBW) configu- 
ration. Multidisciplinary design optimization (MDO) studies of truss-braced wing airplanes suggest 
that optimal designs can have very flexible wings. [2] It is well known that vehicles with long flexible 
wings can have aeroelastic issues related to flutter, gust loads, maneuver loads, limit cycle oscillation, 
ride quality and buckling. [3] The TBW may also have stability issues with low sweep angles and very 
low structural frequencies that can couple with aircraft rigid body modes and flight control systems. 

Other aircraft development efforts that have resulted in very flexible wings are aimed at high en- 
durance. Experience in aircraft design programs intended to improve aerodynamic efficiency such as 
the FliLDA (High Lift-to-Drag Active) wing [4], high altitude long endurance (HALE) [5], and the 
Aeroenvironment Helios crash investigation [6] indicate that gust response requires more thorough 
analysis and validation using nonlinear multidiscipline aeroservoelastic codes. For these reasons it is 
likely that many future projects will necessitate a new level of analysis not seen in current production 
a ire raft design practices. Analyses will include a fluid/structural model capable of simulating poten- 
tially large, and therefore possibly nonlinear deflections. Critical analyses will include flutter and gust 
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loads analyses, while the final design will likely include closed loop flutter suppression and gust load 
alleviation [7,8]. 

Gust analyses have been a standard part of vehicle loads analysis for many decades. Analysis 
to date indicates that gust loads and closed loop gust response of the TBW may be a critical design 
factor. Production vehicle design practice uses gust analysis with linear aerodynamics. Very sophisti- 
cated but fully linear gust models have been developed for both time and frequency domain analyses. 
The Laplace transform of an arbitrary gust wave form has facilitated the development of frequency 
domain models. Alternately, reduced order models developed using time history data of an appro- 
priate parameter set for an aeroelastic model can be imported into a linearized state-space model for 
closed loop analysis. However, the reduced order gust model and the aerodynamic response data 
needed to construct that model must be generated first. Time domain gust analysis has historically 
been performed using a panel code in which the introduction of a perturbation velocity as a local 
angle of attack increment is relatively straight forward. A linear panel code is acceptable if the flow 
field and unsteady aerodynamic response are entirely linear. Linear gust aerodynamics may or may 
not be adequate as more flexible vehicles are designed to fly in the transonic flight regime. If the 
steady state or unsteady flow are nonlinear, a higher order CFD simulation of the gust response will 
be required. 

For this reason there has been an interest in developing gust modeling capability in high fidelity 
nonlinear CFD codes. A technique called the field velocity method (FVM) has been developed and 
is widely being implemented for the simulation of a gust within CFD codes. [9-12]. The present 
paper describes the implementation, verification and use of this gust simulation method within the 
NASA CFD code FUN3D. The first sections briefly describe the FUN3D code and the theoretical 
background for the FVM gust model. These sections arc followed by verification of the FUN3D gust 
model by comparing with theoretical results and with other CFD results. The final sections describe 
a method of creating a reduced order model (ROM) of a gust, and use in simulating a novel sequence 
of gust profiles. 

2 FUN3D Solver 

The Navier-Stokes code FUN3D (fully unstructured Navier-Stokes three-dimensional) is a finite- 
volume unstructured CFD code for either compressible or incompressible flows [13, 14]. Flow vari- 
ables are stored at the vertices of the grid. FUN3D can solve the discrete compressible Euler or 
Reynolds-averaged Navier-Stokes (RANS) flow equations either tightly or loosely coupled with a 
turbulence model on mixed element grids, including tetrahedra, prisms, pyramids and hexahedra. 

The present study uses both the Euler and RANS solution capabilities of FUN3D. The RANS 
simulations include turbulence modeling, performed by loosely coupling the Spalart-Allmaras turbu- 
lence model [15]. Solutions in this study are on either all prismatic or all hexahedral grids. Steady 
state and subiterative solutions are accelerated to convergence by the use of local time stepping [13]. 
Domain decomposition is used to enable distributed parallel computing. 

3 Field Velocity Method of Modeling a Gust 

The FVM physically introduces gusts into a CFD code by utilizing grid velocity [9-12]. Normally 
grid velocity would be associated with physically moving the grid. However, it is possible to introduce 
an arbitrary perturbation velocity in a stationary grid by prescribing the grid velocity at either all or 
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some of the field grid points without actually moving the grid. For instance, the gust profile can be 
defined by a streamwise variation in a perturbation velocity 


for 

and 

for 


u g *(x*,t*) =f x *{t*-x*/Uj) 

v*(x\t*)=f y *{t*-x*/Uj) , 

w g *(x*,t*)=f z *(t*-x*/Uj) , 

f>x*/Uj , 

u g *(x*,t*) = 0 , v g *(x*,t*) = 0 , w g *(x*,t*) = 0 

t* < X* /Uco* . 


( 1 ) 

(2) 

(3) 

(4) 

(5) 

( 6 ) 


The gust profile would be introduced into the flow field by modified nondimensional grid speed 
metrics 

x z i+yT:j + ZT^ = {x T -u g )i+(y T -v g )j + {zz — w g )k . (7) 

The dimensional and nondimensional velocities are related by /aj" where = (u g ,v gl w g ). 

To date this method has been used for relatively simple to moderately complex configurations 
[16]. Further use and development of the method have been performed by Singh and Baeder [11], 
Zaide and Raveh [10, 17], Raveh [16, 18] and Wang et al. [19]. It was recently used to simulate the 
aeroelastic response of a launch vehicle to a sequence of one-minus-cosine gusts [20]. 

This approach has been implemented in the FUN3D CFD code. Gust profiles such as sharp edged 
gusts or one-minus-cosine, as in Figure 1, or even more complex shaped gusts can be introduced [9]. 


4 Discrete Gust Models 


Several discrete gust profiles can be defined in the FUN3D code. They are the Gaussian, the one- 
minus-cosine, the sine and cosine profiles as well as an arbitrary combination of these profiles. Each 
of these gust profiles will be defined here considering the ^-component of gust velocity only, although 
in general all three components of velocity can be defined similarly. 

The first gust to be defined is the Gaussian profile. In nondimensional units, a z-direction gust 
perturbation velocity can be written 


where 


and 


w 


,(x,t) = w, 


- ce 2 


exp*- 


e = 


TMoo 

Lg exp 


c = ln( 2) 


T — 1 tre fexp 


(x-xp) 


( 8 ) 

(9) 

( 10 ) 


The appearance of Mach number in these nondimensional gust parameters is due to the fact that 
FUN3D nondimensionalizes velocities using free stream speed of sound whereas the gust itself con- 
vects at the free stream velocity. Figure 2 illustrates the Gaussian gust profile. The requirement that 
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the profile have the width shown in Figure 2 makes the definition of the constant c in equation 9 
clear. Table 1 relates these parameters to FUN3D namelist input. Note the minor nomenclature in- 
consistency, that the Gaussian profile is given a subscript ”exp” in this discussion and in the FUN3D 
namelist input. 

A full cycle one-minus-cosine profile is frequently used in discrete gust analyses. In the present 
formulation the term one-minus-cosine will be used to denote a half cycle one-minus-cosine followed 
by a hold function. Elsewhere this has been called a ramp-hold function. In terms of the present 
definition, a full cycle one-minus-cosine profile can be constructed by combining two of the current 
profiles with equal magnitudes but opposite signs separated by a half cycle. The present definition also 
allows for the combination of multiple one-minus-cosine functions into a smoothly varying profile. 
Using the present definition, in nondimensional units, the equation for a one-minus-cosine gust profile 
z-component of velocity is 

Wg(x,t) = -w omc [l -cosG] (11) 

where 

Q = nzM^ 0 <t<%- (12) 

L Somc M °° 

and 

6 = % for T > . (13) 

The nondimensional time parameter is 

T = r-r„, , (14) 

rt J omc v ' 

Figure 3 illustrates the one-minus-cosine gust profile, while Table 1 relates these parameters to 
FUN3D namelist input. 

In nondimensional units, the equation for the sine gust profile is 


Wg(x,t) = WsinSind 


(15) 


where 


and 


2kxM„ 

Lg sin 


for X > 0 


T — f tre f sin 


(x-xo) 

Moo 


(16) 

(17) 


Figure 4 illustrates the sine gust profile, while Table 1 relates these parameters to FUN3D namelist 
input. 

For the cosine gust profile, one has 


w g (x,t) = w cos cos6 


(18) 


The cosine gust profile parameters 0 and x are defined in the same way as for the sine gust profile. 

To simulate gusts, the only requirement is that the FUN3D code must be operated in unsteady or 
time accurate mode. As shown in Figure 1 , a forward starting point Go) for the gust velocity can be 
specified. The default is xo = 0. Figure 5 shows a typical FUN3D gust data namelist input. Note that 
as defined in each of these four profiles and in equations 1-6, gust shapes vary in the v-dircction only. 
The gust velocity is uniform in the y and z directions. 


5 Verification 


The functioning of the FVM is verified by comparison with other published results, and with theo- 
retical results. Parameswaran and Baeder compute the response of an NACA 0006 airfoil to indicial 
angle of attack change [12]. The approach incorporates the step change in angle of attack as a step 
change in grid velocity over the entire flow domain instead of only at the airfoil surface. Note that an 
angle of attack change is different than a gust in that a gust includes a down stream convection. For 
the present simulation the instantaneous angle of attack change was incorporated by setting a uniform 
change in grid velocity throughout the computational domain. Even though it is not strictly a gust, it 
serves as a valuable test because there are theoretical results for a step angle of attack change against 
which the current model can be checked. The data of Parameswaran and Baeder is generated using 
the TURNS code [21]. The finest mesh they used was a C-type structured inviscid mesh with 251 
grid points in the chordwise and 61 grid points in the direction normal to the airfoil surface. 

The present grid was created with the AFLR3 grid generator [22]. It is an inviscid two-dimensional 
grid with 52,000 cells. Figure 6 shows the time history of normal force coefficient (C^/Aa) due to 
a step angle of attack, Aa. The non-dimensional time parameter S is defined S = 2UJ*t* /C * . The 
figure presents the responses at Mach numbers 0.3 and 0.5. The force coefficient response time his- 
tories presently computed with FUN3D compare very well at both Mach numbers with the responses 
computed by Parameswaran and Baeder. 

Figure 7 shows the normal force coefficient per angle of attack response for short time. The 
present results are compared with the computed results of Parameswaran and Baeder and the ana- 
lytical result of Lomax [23]. The expression for normal force coefficient at small times is given by 
Lomax as 


C N /\a = {-^-){\ 


(1 - Moo) 
2 Mec 


S} 


(19) 


for 


0 < S < 2Afoo/(l +MJ) 


(20) 


The CLD results show oscillations in the initial few time steps. These oscillations may be numerical 
noise generated by the introduction of a discontinuity in grid velocity at the start of the simulation. 
Disregarding the early oscillations, the CLD response remains linear during the time period over 
which the theory says it should be linear. This comparison is a good confirmation of the present 
implementation of the grid speed metrics. 

The next verification tests compare LUN3D responses to a 5 chord length and a 25 chord length 
full cycle one-minus-cosine gust profile and a 25 chord length sine gust profile with that computed by 
Zaide and Raveh [10]. This and the following cases are true gust simulations as defined in equations 
1-6. The configuration is a NACA 0012 two-dimensional airfoil. The condition is at Mach 0.2 and 
a = 0 degrees. The simulations of Zaide and Raveh were done with a C-type inviscid mesh having 
399 grid points in the chordwise direction and 7 1 grid points normal to the wing surface. In the present 
computations, an unstructured prismatic grid was constructed using the grid generator ALLR3. It had 
52,000 grid points. Inviscid simulations were performed. Three gust calculations as defined above 
were performed. In each example the gust velocity has z-component only. The gust velocity function 
magnitudes and gust lengths are as defined in Zaide and Raveh [10]. Ligure 8 shows the 5 chord length 
gust Cl response while Ligure 9 shows the response to a 25 chord length gust. Despite differences 
in the meshes and solution methodologies, the present results and those of Zaide and Raveh compare 
well. 


9 


Finally, the airfoil is subject to a single cycle of a sine function profile. The magnitude and fre- 
quency of the sine function arc defined in Zaide and Raveh [10]. Beyond the first cycle of oscillation 
the gust velocity is zero. The lift coefficient response to a sine wave gust profile is shown in Figure 
10. Once again, the comparison of the present FUN3D lift coefficient response time history with that 
of Zaide and Raveh is very good. In all three cases, despite the differences in grid type, the response 
computed with FUN3D compares reasonably well with the response computed by Zaide and Raveh. 


6 Development of an ARMA Reduced Order Gust Model 

Having verified that the FUN3D gust model performs as expected for several test cases, its use will 
be illustrated by several examples. Here, the focus will be on extracting reduced order gust models. 
Reduced order models of gusts generated from CFD gust system identification allow a variety of 
novel, or previously unseen, sequences of gust profiles to be simulated in an efficient open or closed 
loop state space model of an aircraft. For example, Raveh used a reduced order gust model to sim- 
ulate gust response of an aircraft [9]. An ARMA model was constructed of the CFD aerodynamic 
response to the gust which was then used in a state space model of the aeroelastic vehicle response. 
More complicated gust solutions were then computed using the Federal Aviation Regulations (FAR) 
required gust definition by power spectral density (PSD) based on the Dryden gust shape filter [9]. 
Other approaches such as a convolution integral have been used to create reduced order models of 
gusts [19]. 

This section outlines a reduced order model of a gust that uses the Auto-Regressive Moving- 
Average (ARMA) method. An ARMA model is intended to recreate an output of interest. To compute 
the output {y} at time step /, the model is written 

p Q 

y\ = Yj a kyi-k + Yj b k w si~k ( 21 ) 

k = 1 k=0 

The output can be e.g. yj = (C+C,„, {G\ T ) where {G\ is a vector of generalized forces. There are a 
variety of methods available to compute the coefficients a and b, such as the fast orthogonal search 
method [24], group method of data handling [25], or more recently the optimal parameter search 
(OPS) method based on affine geometry developed by Lu et al. [26]. In the present paper, the arrays 
of coefficients a and b will be obtained by the OPS ARMA method. In that method, the parameters 
P and <2+1 represent the maximum AR and MA model orders, respectively. The actual number of 
terms is usually less than these orders. The value of this method is that it obtains the ARMA model 
with the minimum or optimal set of terms. 

To illustrate the use of this method, a reduced order model of the response of Cl and Cm to a gust 
for the NACA 0012 airfoil is developed. The condition is at Mach 0.65, a = 2 degrees. The Reynolds 
number is 1 . 1 million based on chord length. In this case a two-dimensional unstructured hexahedral 
viscous grid was created from a structured C-type mesh with 305 grid points in the chordwise direc- 
tion and 129 points normal to the wing surface. For the excitation of the flowfield, a Gaussian doublet 
pulse, shown in Figure 1 1 , was used. The Gaussian doublet pulse is composed of two Gaussian pro- 
files, the second having a slight time lag compared to the first. The gust excitation induces responses 
in lift and moment coefficients. Time histories of the gust response component of the lift and moment 
coefficients are shown in Figure 12. From those responses, an ARMA model with a maximum of 5 
AR and 5 MA terms was created. 


10 


To test the validity of the ARMA model, a novel gust input, previously unseen by the ARMA 
model, is input into the system. The novel gust input is shown in Figure 13. This gust input is 
a sequence of 20 one-minus-cosine profiles of varying lengths and amplitudes. At the end of the 
sequence the gust amplitude returns to zero. Figure 14 shows the direct FUN3D simulated gust 
responses and the ARMA model response to this new input. The excellent agreement between the 
FUN3D and ARMA model lift and moment coefficients validates the approach. Note that the lift and 
moment coefficient error sizes are similar; the scale of the plots makes the moment coefficient error 
appear to be larger. 


7 Aeroelastic Simulation Using a POD/ARMA Reduced Order Gust 
Model 


Rather than creating an ARMA model capable of predicting only integrated coefficient time histories, 
it would be of interest to create a reduced order model that is capable of predicting the time history 
response of the distributed surface pressures to a gust input. Such a model can be constructed using a 
POD of the covariance matrix of the unsteady pressure coefficients. In this model, the POD eigenvec- 
tors provide the spatial distribution of the principle components of the unsteady pressure response. 
An ARMA model can be developed that will provide the time varying component of the model. 
The combined POD eigenvectors and ARMA model coefficients provide a method of predicting the 
unsteady surface pressure responses to an arbitrary gust. 

A Gaussian doublet (zero mean) gust velocity profile is used to generate the pressure response 
training data. The covariance matrix of the dynamic component of surface pressure coefficients is 
defined 

R = [A] [X] r (22) 

where 

[*f=[Si ••• S, ••• %]. (23) 

is a matrix the columns of which are composed of N f snap shots of the dynamic component of the 
surface pressure coefficient distribution. The snap shot at the step I is 




Ac 


Pl,l 


Ac 


PNJ. 


(24) 


The N t snap shots are the training data for the following POD and ARMA models. The POD modes 
are the eigenvectors associated with the eigenvalue problem X T X<$ = <E>A where <f> is the matrix of 
eigenvectors and A are the eigenvalues. The size of the matrix <P is reduced to form a truncated POD- 
basis <I>, with dimension N s x Npod by retaining only the eigenvectors associated with the largest 
eigenvalues. 

The training data can be written 

S/ = [<U] A ■ (25) 

The term /3/ is an Npod dimension array of coefficients at time step /.The coefficients /3/ can be found 
by several methods. One method is by finding the solution /3 that gives the minimum 7/ 


Ji= min |‘4> r j3/-E/| 2 

fieR N POD 


(26) 
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Alternately, since the covariance matrix R is symmetric, and if the eigenvector matrix <t>, is scaled to 
have a unit norm, then 

ft = [<ftf 2, . (27) 

The former approach is taken here. 

Corresponding to equation 25, the pressure coefficient distribution at time step / predicted by the 
model to be developed is 

S/ = [<ft]ft (28) 

where the predicted dynamic component of the surface pressure coefficients is 




Acpi,/ 

_^CpN s ,l 


(29) 


and ft is a predicted array of coefficients at time step 1. A predictive ARMA model can be created to 
compute the array sequence [ft , .... fty, ]. At time step / 


p „ Q 
ft = L N k$l-k + L b kWg,- k 

k= 1 k = 0 


(30) 


where the diagonal matrices a k and arrays b k contain the autoregressive and moving average terms of 
the ARMA model. The OPS ARMA method, discussed in the last section, provides the estimate of 
the coefficients a and b. Equations 30 and 28 provide the means to predict the pressure distribution 
due to a novel gust input 

p , Q 

%=Y [‘ft] [ a ]kPl-k + Y, [‘ft ] b k W gl-k ■ ( 31 ) 

k= 1 k = 0 

This model reduces the N t training snap shots of N s surface pressure coefficients to a predictive model 
having (N s +P + Q + 1 ) x Npod terms. The predicted surface pressure distribution can be integrated 
to provide the generalized force input due to gust. The generalized force due to gust can be written 


G g = q4<l>} T [s}{Y [^ MkPl-k+Y i®r]b k w gl _ k } 


k= 1 


*= 0 


(32) 


where [0] is a 3 N s x N m matrix 


0v. 1 1 

' ftdA„ 

ft, 11 • 

‘ ft,lA m 

ft, 11 • 

' ftdA„ 

ft Asl ■ 

• ftcA *N m 

ftdYd ' 

• ftAAm 

.ft Ad ' 

• ftAAm 


(33) 


The integration matrix 5 with dimension 3 N s x N s is defined so that 


f = q 00 [s}£ 


(34) 
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where the nodal force array / is given by 


' fx, l" 

fy, l 

fz, 1 

fx, N s 

fy, N s 

Jz,N s . 


( 35 ) 


The combined aeroelastic and gust models are cast in state space form. The unsteady aerody- 
namic terms due to modal displacements in the aeroelastic model are constructed using the Roger 
approximation [27]. The state space model with lag states % and gust input G „ , is given by 


m=[A]M+G, , 


(36) 


where 


[A] 


0 / 0 

—M~ X Q. -AT 1 A M~ 1 q 00 d . 
0 e 7 


(37) 


The Roger approximation is used to create the e, 7 matrices and the d coefficients. The Roger model 
is created by exciting each of the modes separately using, in this case, Gaussian pulses. With the 
complete database assembled, the numerical integration of equation 36 can be performed. The present 
simulations are computed with a 5 th order implicit Euler backward difference scheme. 

As an example illustrating the use of this method, aeroelastic simulations of the Benchmark Active 
Controls Technology (BACT) aeroelastic model will be performed with and without gust input. The 
BACT aeroelastic model is composed of a rigid unswept wing with constant chord along the span. 
The chordwise section is an NACA 0012 airfoil. In the present CFD model this wing is approximated 
as a two-dimensional airfoil. Thus, three-dimensional wing tip effects are absent in the present CFD 
model. In the wind tunnel BACT model, the wing root is attached to an apparatus that allows elastic 
plunge and pitch degrees of freedom. The elastic and mass properties of the model are found in 
Table 2. Based on this model the plunge and pitch modes are constructed. The two modes are shown 
in Figure 15. The vertical modal amplitudes are shown. Note that the first mode is predominantly 
plunge, while the second mode is predominantly pitch. These modes are used in the CFD aeroelastic 
model and in the integration of the gust surface pressures. After creating the Roger approximation and 
gust model, a state space formulation of the aeroelastic model can be developed. The gust forcing 
to the state space model will be in the form of generalized force integrated from surface pressure 
coefficients which in turn have been generated by the gust input. 

The BACT case to be simulated is case SEFC2 of reference [28]. This case is at Mach 0.747, 
a = 1 .78 degrees in R-12 heavy gas. The Reynolds number is 1.1 million based on chord length. The 
experimental flutter onset is at a dynamic pressure of 15 1 .6 pounds per square foot (psf). Six lag states 
are used to simulate this case. The aeroelastic response is initiated by setting initial modal velocities 
for each mode. To test the accuracy of the Roger approximation, the dynamic pressure at which 
flutter onset occurs was simulated first with no gust input. Flutter onset is taken to be the dynamic 
pressure at which the oscillatory amplitude of either the pitch or plunge mode is constant. Above that 
dynamic pressure, the oscillatory amplitude grows. Flutter onset was approached by step increases 
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in the dynamic pressure by several psf per solution, starting from q „ = 60 psf. The calculated flutter 
onset (with G g = 0) is approximately 132 psf. While about 15 percent below the experimental flutter 
onset, this value of dynamic pressure is sufficiently close for the puipose of the present demonstration. 

The time histories of the two BACT model modes (with G g = 0) at a dynamic pressure of 60 
psf are shown in Figure 1 6. The symbol (A) indicates that it is the dynamic component only. Two 
solutions are shown in that figure. One solution is a time accurate aeroelastic solution using FUN3D. 
The other solution is generated using the state space model. With the exception of slight over shoots 
in the state space solution early in the simulation, the time histories of the generalized displacements 
computed using FUN3D and using the state space model are in close agreement. 

A simulation is next performed that includes a gust input. A sequence of 7 sine functions in z- 
velocity (w g ) are used. The gust velocity time history is shown in Figure 17(a). This sequence of sine 
functions is used to create a pseudo-random vertical continuous gust environment. The amplitudes of 
the sine functions are set to roughly approximate a 0.5 — 1.0 percent free stream turbulence level. The 
generalized forcing due to this gust input is integrated from the pressures computed using equations 
28- 30. The 12 leading POD eigenvectors are used to form the truncated basis <£>,.. The gust ARMA 
model matrix a and array b each have dimension 8. The simulations are again initiated with modal 
velocities. The resulting time histories of the dynamic component of generalized displacements are 
shown in Figure 17(b) and 17(c). The FUN3D simulated time history and the state space with 
ARMA gust input time history agree very well. Figure 18 shows the leading and trailing edge 
dynamic displacements for the case with gust, computed with the state space model. With respect to 
a 16 inch wing chord, the vertical displacements shown are very small. 


8 Concluding Remarks 


This paper has presented the implementation of gust modeling capability using the field velocity 
method into the FUN3D CFD code. A detailed discussion is presented of the various gust profiles 
available and the corresponding FUN3D input. This paper serves as a users manual for FUN3D gust 
analyses using a variety of gust profiles. The gust capability has been verified by computing the 
response of an airfoil to a step angle of attack change and comparing this result with theoretical and 
other computational results. The present simulations have been shown to agree very well with the 
theoretical and other CFD gust simulations. Development of an Auto-Regressive Moving-Average 
(ARMA) reduced order gust model is also presented. The system identification required to develop 
the ARMA model is performed with a Gaussian gust profile in the FUN3D code. To verifiy the 
utility of an ARMA reduced order model approach, the response of the system to novel gust inputs 
is shown. ARMA simulated response to a sequence of one-minus-cosine gusts is shown to compare 
well with the same gust profile computed with FUN3D. Proper Orthogonal Decomposition (POD) 
is combined with the ARMA modeling technique to provide a means of obtaining the time varying 
pressure coefficient increment distribution due to a novel gust profile. This approach is combined 
with an aeroelastic reduced order model. Simulation of the aeroelastic response of a pitch/plunge 
airfoil to a gust environment is computed with this reduced order model. This result is compared with 
a direct simulation of the system in the FUN3D code. The two results are found to agree very well. 
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Table 1: Gust parameter and FUN3D namelist ” gust data” definitions. 


Namelist 

Name 

Para- 

meter 

Gust Profile 

Definition 

l_gust_cos 

Lgcos 

Cosine 

Wave length of cosine profile 

l_gust_exp 

Lg exp 

Gaussian 

Half-height length of Gaussian profile 

l_gust_omc 

Lgomc 

One-minus-cosine 

Length of one-minus-cosine profile 

l_gust_sin 

Lg sin 

Sine 

Wave length of sine profile 

ngust_cos 

Vgcos 

Cosine 

Number of cosine profiles 

ngust_exp 

Ng exp 

Gaussian 

Number of Gaussian profiles 

ngust_omc 

N 

y gome 

One-minus-cosine 

Number of one-minus-cosine profiles 

ngust_sin 

Ng sin 

Sine 

Number of sine profiles 

tref _gust_cos 

tre fcos 

Cosine 

Reference time for cosine profile 

tref _gust_exp 

^ re f exp 

Gaussian 

Reference time for Gaussian profile 

tref _gust_omc 

tre fomc 

One-minus-cosine 

Reference time for one-minus-cosine profile 

tref _gust_sin 

tre f sin 

Sine 

Reference time for sine profile 

u_gust_cos 

Mcos 

Cosine 

x velocity magnitude of cosine profile 

u_gust_exp 

Mexp 

Gaussian 

x velocity magnitude of Gaussian profile 

u_gust_omc 

Momc 

One-minus-cosine 

x velocity magnitude of one-minus-cosine profile 

u_gust_sin 

Usin 

Sine 

x velocity magnitude of sine profile 

v_gust_cos 

V cos 

Cosine 

y velocity magnitude of cosine profile 

v_gust_exp 

V exp 

Gaussian 

y velocity magnitude of Gaussian profile 

v_gust_omc 

Vomc 

One-minus-cosine 

y velocity magnitude of one-nrinus-cosine profile 

v_gust_sin 

V sin 

Sine 

y velocity magnitude of sine profile 

w_gust_cos 

Wcos 

Cosine 

z velocity magnitude of cosine profile 

w_gust_exp 

Wexp 

Gaussian 

z velocity magnitude of Gaussian profile 

w_gust_omc 

W omc 

One-minus-cosine 

z velocity magnitude of one-nrinus-cosine profile 

w_gust_sin 

W 'sin 

Sine 

z velocity magnitude of sine profile 

xO_gust 

x 0 

All 

x starting location of gust perturbation velocity 
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Table 2: BACT model parameters. 


Parameter 

Name 

Dimensions 

Value 

M 

Mass 

slugs 

5.966 

S a 

Static Offset 

slug- ft 

0.0142 

la 

Pitch Moment of Inertia 

slug- ft 2 

2.8017 

K h 

Plunge Stiffness 

Ib/ft 

2659 

K a 

Pitch Stiffness 

lb — ft / rad 

2897 
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&gust_data 

ngust_sin = 2 
ngust_cos = 1 
ngust_omc = 1 
xO_gust = 0. 
w_gust_sin(1) =0.0005 
w_gust_sin(2) =-0.0005 
w_gust_cos(1 ) = 0.0100 
w_gust_omc(1)= 0.0010 
l_gust_sin(1) =2500. 
I_gust_sin(2) = 2000. 
l_gust_cos(1) = 1500. 
l_gust_omc(1) = 1000. 
tref gust sin(l) =0. 
tref_gust_sin(2) = 0. 
tref_gust_cos(1 ) = 0. 
tref_gust_omc(1 )= 0. 

/ 


Figure 5: FUN3D namelist gust input. 
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Figure 6: Cn response (per Aa) to a step change in angle of attack. 
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Figure 7: Cn response (per Aa) to a step change in angle of attack. 
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Figure 8: Responses due to one-minus-cosine gust, 5 chord gust length. 



Figure 9: Responses due to one-minus-cosine gust, 25 chord gust length. 
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Figure 10: Responses due to a single-period sine gust, 25 chord gust length. 




0.4 


FUN3D time marching 




(b) Pitch moment A Cm time history. 


Figure 12: Force and moment coefficient responses to Gaussian doublet gust input, Mach 0.65, a = 2 
degrees. 
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(a) Force coefficient A Q. time history. 



(b) Pitch moment A Cm time history. 


Figure 14: Force and moment coefficient responses to a sequence of 20 one-minus-cosine gust pro- 
files, Mach 0.65, a = 2 degrees. 
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(b) Mode 2 generalized displacement time history. 


Figure 16: BACT dynamic component of generalized displacements, no gust, Mach 0.747, q <*, = 60 
psf, a = 1.78 degrees, R-12 heavy gas. 
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(a) Gust velocity w g time history. 
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(b) Mode 1 generalized displacement time history. 



t ~ sec. 

(c) Mode 2 generalized displacement time history. 
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Figure 17: BACT dynamic component of generalized displacements, with gust, Mach 0.747, = 60 

psf, a = 1.78 degrees, R-12 heavy gas. 




(a) Leading edge displacement time history. 



(b) Trailing edge displacement time history. 

Figure 18: BACT dynamic component of z displacements due to gust computed with state space 
model. Mach 0.747, q„ = 60 psf, a = 1.78 degrees, R-12 heavy gas. 
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